% 有限差分求解平面弹性力学问题...
% 目前只支持设置位移边界
% 可能存在bug！！！仅供娱乐

clc
clear

L = 5;
dx = 0.1;

[x y] = meshgrid(0:dx:L);
n = size(x, 1);


E = 10;
nu = 0.3;
G = E/(2*(1+nu))
lambda = nu*E/((1-2*nu)*(1+nu))

%A = zeros(n,n);
A = spalloc(2*n^2, 2*n^2, 10*n^2); %根据AI建议，使用稀疏矩阵优化内存
b = zeros(2*n^2,1);

boundary_right = @(x,y) 0.1*sin(pi/L*y).^10; %右边界

printf('Setting up Stiffness Matrix...\n')
for i = 1:n
    for j = 1:n
        _k = i+(j-1)*n;

        if i == 1 || i == n || j == 1 || j == n %如果是边界
            A(_k,_k)=1;
            b(_k) = 0;

			if j == n %右u1位移边界
                b(_k) = boundary_right(x(i,j),y(i,j));
            end
        else
            A(_k,_k) = -4*G;
            A(_k,_k + 1) = G;
            A(_k,_k - 1) = G;
            A(_k,_k + n) = G;
            A(_k,_k - n) = G;

			A(_k, _k) += -2*(G+lambda);
			A(_k, _k + n) += (G+lambda);
			A(_k, _k - n) += (G+lambda);

            A(_k,_k+n^2+n+1) += 1/4*(G+lambda);
            A(_k,_k+n^2-n-1) += 1/4*(G+lambda);
            A(_k,_k+n^2-n+1) += -1/4*(G+lambda);
            A(_k,_k+n^2+n-1) += -1/4*(G+lambda);

        end
    end
end

for i = 1:n
    for j = 1:n
        _k = n^2+i+(j-1)*n;

        if i == 1 || i == n || j == 1 || j == n %如果是边界
            A(_k,_k)=1;
            b(_k) = 0;
        else
            A(_k,_k) = -4*G;
            A(_k,_k + 1) = G;
            A(_k,_k - 1) = G;
            A(_k,_k + n) = G;
            A(_k,_k - n) = G;

			A(_k, _k) += -2*(G+lambda);
			A(_k, _k + 1) += (G+lambda);
			A(_k, _k - 1) += (G+lambda);

            A(_k,_k-n^2+n+1) += 1/4*(G+lambda);
            A(_k,_k-n^2-n-1) += 1/4*(G+lambda);
            A(_k,_k-n^2-n+1) += -1/4*(G+lambda);
            A(_k,_k-n^2+n-1) += -1/4*(G+lambda);
        end
    end
end

_u = zeros(2*n^2,1);
%_u = inv(A)*b;
printf('Solving Equations...\n')
_u = gmres(A,b,20,1e-4,100); %使用迭代法求解会比求逆快得多

_u1 = _u(1:n^2);
_u2 = _u(n^2+1:2*n^2);
u1 = reshape(_u1,n,n); %将u变回n*n的场
u2 = reshape(_u2,n,n); %将u变回n*n的场

figure();
hold on
axis equal
quiver(x,y,u1,u2)
xlabel('x');
ylabel('y');

